void vscale(caja *bx)
{
	///Termostato empleando el "Momentum Scalling Method"

	///Variables externas: Ahi que ver de donde vienen y como generarlas
	double atemp = bx->temp;
	double temp=0.0;
	double taut=5000.0;
	////Variables de entrada del termostato
	double alpha;
	for(int j=0;j<3;j++)
	{
		for(int i=0;i<bx->npart;i++)
		{
			temp+= (bx->part[i].vel[j])*(bx->part[i].vel[j]);
		}
	}
	temp*=1.0/(3.0*bx->npart*kB);
	//alpha = sqrt(1.0+(atemp/temp-1.0)/taut);
	alpha = sqrt(atemp/temp);

	//alpha = sqrt(atemp/temp);
	//alpha=1.0;
	cout<<"alpha: "<<alpha<<endl;
	cout<<"atemp: "<<atemp<<"  "<<temp<<endl;
	bx->alpha= alpha;
}

void vscale2(caja *bx)
{
	///Termostato empleando el "Momentum Scalling Method"

	///Variables externas: Ahi que ver de donde vienen y como generarlas
	double atemp = bx->temp0;
	double temp=0.0;
	////Variables de entrada del termostato
	double alpha;
	for(int j=0;j<3;j++)
	{
		for(int i=0;i<bx->npart;i++)
		{
			temp+= (bx->part[i].vel[j])*(bx->part[i].vel[j]);
		}
	}
	//temp*=1.0/(3.0*bx->npart*kB);
	temp*=1.0/(3.0*bx->npart);
	alpha = sqrt(atemp/temp);
	for(int i=0;i<bx->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			bx->part[i].vel[j]*=alpha;
		}
	}
	//cout<<"atemp: "<<atemp<<"  "<<temp<<endl;
}

void dscale(caja *bx)
{
//	///Manostato modificando el volumen de simulacion y reescalando posicion de
//	///las particulas
//
//	///Variables externas: Ahi que ver de donde vienen y como generarlas
//	double beta, delta, apres,pres,taup, comp, facpot;
//	//////////////////
//      int nfi,nacc;
//      beta=pow(1.0-comp*delta*(apres-pres)/taup,1.0/3.0);
//
//      for(int i=0;i<3;i++)
//      {
//    	  bx->lbox[i]*=beta;
//    	  for(int j=0;j<bx->npart;j++)
//    	  {
//    		  bx->part[j].pos[i]*=beta;
//    	  }
//      }
//      bx->vol=bx->lbox[0]*bx->lbox[1]*bx->lbox[2];
}


